%plots perturbation population response
%plots population resonse on running onset (fb+dark)

th_run=0.05;
plotit=0;
fb_ac=[];pb_ac=[];dark_r_ac=[];dark_b_ac=[];
for site=1:size(proj_meta,2)
    no_time_points = size(proj_meta(site).rd,2);
    
    act=[];
    type_sessions={'f','p','d'};
    for tp=1:no_time_points
        
        act=[];
        for lyr=1:size(proj_meta(site).rd,1)
            act=[act; proj_meta(site).rd(lyr).act];
        end
        avg_act=mean(act);
        
        sess_start=proj_meta(site).rd(lyr).nbr_frames;
        sess_start=cumsum(sess_start);
        sess_start=[1 sess_start+1];
        ac_sess=find(strcmp(proj_meta(site).rd(lyr).session,type_sessions(1))==1);
        indices=[];
        for kk=1:length(ac_sess)
            indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
        end
        
        velM = proj_meta(site).rd(lyr).velM_smoothed;
        velP = proj_meta(site).rd(lyr).velP_smoothed;
        visSpeed = proj_meta(site).rd(lyr).vis_speed;
        run=velM/max(velM)-mean(velM);
        vis=velP/max(velP)-mean(velP);
        velM_r = find(run>th_run);
        velP_r = find(vis>th_run);
        if plotit
            figure;
            hold on
            plot(avg_act,'k'), hold on
            plot(run+1.5,'r')
            plot(vis+2.6,'g')
        end
        
        indices = intersect(indices,velM_r);
        
        fb_ac(site) = mean(avg_act(indices));
        if plotit, plot(indices,ones(1,length(indices))*.9,'r.'),end
        run1(site) = mean(run(indices));
        
        ac_sess=find(strcmp(proj_meta(site).rd(lyr).session,type_sessions(2))==1);
        indices=[];
        for kk=1:length(ac_sess)
            indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
        end
        ind1 = intersect(indices,velP_r);
        pb_ac(site) = mean(avg_act(ind1));
        run2(site)=mean(run(ind1));
        run2b(site)=mean(run(intersect(indices,velM_r)));
        if plotit, plot(ind1,ones(1,length(ind1))*.91,'g.'),end
        
        ac_sess=find(strcmp(proj_meta(site).rd(lyr).session,type_sessions(3))==1);
        indices=[];
        for kk=1:length(ac_sess)
            indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
        end
        ind1=intersect(indices,velM_r);
        dark_r_ac(site) = mean(avg_act(ind1));
        if plotit, plot(ind1,ones(1,length(ind1))*.89,'b.'),end
        run3(site)=mean(run(ind1));
        ind1=intersect(indices,find(run<0.005));
        dark_b_ac(site) = mean(avg_act(ind1));
        if plotit, plot(ind1,ones(1,length(ind1))*.88,'k.'),end
        title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
    end
end

% average act- across all animals
result = ([mean(fb_ac) mean(pb_ac) mean(dark_r_ac) mean(dark_b_ac)]-1)*100;
sem_r = ([SEM(fb_ac) SEM(pb_ac) SEM(dark_r_ac) SEM(dark_b_ac)])*100;
figure,errorbar(1:4,result,sem_r,'o');
ylim([0 11])
xlim([0 5])
set(gca,'XTick',(1:4))

result = [mean(run1) mean(run2) mean(run3)];
sem_r = [SEM(run1) SEM(run2) SEM(run3)];
figure,errorbar(1:3,result,sem_r,'o');
hold on 
errorbar(2,mean(run2b),SEM(run2b),'o');
set(gca,'XTick',(1:3))
xlim([0 5])
ylabel('average running (AU)')

